Tracking the precession of compact binaries from their gravitational- wave signal 
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We present a simple method to track the precession of a black-hole-binary system, using only information 
from the gravitational-wave (GW) signal. Our method consists of locating the frame from which the magnitude 
of the (£ = 2, \m\ = 2) modes is maximized, which we denote the "quadrupole-aligned" frame. We demonstrate 
the efficacy of this method when applied to waveforms from numerical simulations. In the test case of an 
equal-mass nonspinning binary, our method locates the direction of the orbital angular momentum to within 
(A0,A<p) = (0.05°, 0.2°). We then apply the method to a q = M 2 /Mx = 3 binary that exhibits significant 
precession. In general a spinning binary's orbital angular momentum L is not orthogonal to the orbital plane. 
Evidence that our method locates the direction of L rather than the normal of the orbital plane is provided by 
comparison with post-Newtonian (PN) results. Also, we observe that it accurately reproduces similar higher- 
mode amplitudes to a comparable non-spinning (and therefore non-precessing) binary, and that the frequency of 
the (ft = 2, \m\ = 2) modes is consistent with the "total frequency" of the binary's motion. The simple form of 
the quadrupole-aligned waveform will be useful in attempts to analytically model the inspiral-merger-ringdown 
(IMR) signal of precessing binaries, and in standardizing the representation of waveforms for studies of accuracy 
and consistency of source modelling efforts, both numerical and analytical. 

PACS numbers: 



I. INTRODUCTION 

Black-hole-binary mergers are expected to be key sources 
for gravitational- wave (GW) astronomy Q]. Accurate theo- 
retical models of the GW signal are necessary to both detect 
these sources and to determine their physical parameters and 
their location in the universe. The GW signal from the inspiral 
can be calculated by analytic approximation techniques J2][3], 
and the merger of the two black holes and ringdown of the 
final black hole can be calculated from numerical simulations 
in full General Relativity HQ. 

Numerical simulations can produce waveforms for only 
discrete points in the parameter space of binary configura- 
tions, but significant progress has been made in synthesiz- 
ing information from post-Newtonian (PN) and effective-one- 
body (EOB) methods, numerical relativity (NR), and pertur- 
bation theory, to produce analytic models of the complete 
inspiral-merger-ringdown signal over some regions of the pa- 
rameter space. Most models to date treat nonspinning bina- 
ries EHTHl, or binaries in which the black-hole spins do not 
precess |[T9l - l2"Tl (although there has been one first attempt at 
a precession model ll22l '). 

Precession adds a number of complications. When the spins 
are not parallel to the orbital angular momentum their orien- 
tation varies with time, as does the orbital angular momentum 
itself; the orbital plane precesses. Both the precession of the 
spins and of the orbital plane each introduce modulations into 
the GW amplitude, oscillations into the GW frequency, and 
variations in the distribution of signal power across different 
harmonics of the waveform. All of these complicate efforts to 
produce an analytic model of precessing-binary waveforms. 
In addition, they make it difficult to uniquely characterize the 
wave signal. For example, the total phase of the dominant 
mode of the signal depends on the initial orientation of the or- 



bital plane. This makes it difficult to determine if two wave- 
forms were produced by the same binary configuration, or to 
compare independent numerical simulations, a task that is rel- 
atively simple for non-precessing non-eccentric binaries ll23l - 

We propose a method to put a precessing-binary waveform 
into a particularly simple form. The method is based on find- 
ing a preferred time dependent coordinate system for the grav- 
itational wave signal, which tracks the precession. 

Gravitational wave signals are most conveniently expressed 
in terms of spherical harmonics of spin-weight s = —2, 
7^(0, <p), where (9,(p) are the standard polar coordinates 
on the unit sphere. The dominant modes are the quadrupole 
modes, where 1 = 2, and —2 < m < 2. If the system is rotated, 
the modes of a particular I mix among each other according 
to the transformation law described in Appendix [A] 

As can be seen from standard post-Newtonian and post- 
Minkowskian descriptions, binary systems emit gravitational 
waves predominantly in the direction orthogonal to the orbital 
plane. Correspondingly, if our system is oriented such that this 
direction is along the z-axis, then we expect that the dominant 
signal is given by the (I = 2, \m\ = 2) spherical harmonics 
of the wave. The modes \m\ = 1 vanish when the two black 
holes can be exchanged by symmetry, and m = is a non- 
oscillating mode related to memory effects, see e.g., ||26ll27l . 
If we choose different (rotated) coordinates (9' , <p') to define 
a new basis Yj L(9', <p'), then mode mixing will complicate the 
spherical harmonic description of the signal, and for example 
even an equal mass nonspinning binary will exhibit nonvan- 
ishing \m\ = 1 modes. We illustrate this effect in Sec. IV A 



Therefore, we can determine a preferred direction from the 
wave signal alone by finding the orientation that maximizes 
the (£ = 2, |m| = 2) modes. This is the method that we will 
discuss in this paper, and we will refer to waveforms that are 
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given in terms of spherical harmonics that are aligned with 
this direction as "quadrupole-aligned" waveforms. 

In a precessing system there are two contributions to the 
frequency of the binary motion: the frequency of the motion 
about the orbital-plane axis, ft) rb> which will increase during 
a non-eccentric inspiral as a monotonic function, and the fre- 
quency of the motion due to the precessional motion, which 
will oscillate as a function of time. The total frequency of the 
motion is co = (O ol ±, — cpcos 9, where 9 is the inclination of the 
normal to the orbital plane from the z-axis, and <p is the rota- 
tion of the normal about the z-axis in the xy plane. (This cor- 
responds to the result in Eq. (3.10) in ll28ll .) In a kinematical 
description of the binary, these two frequencies together pre- 
scribe the bodies' acceleration, which is the dominant source 
of gravitational radiation. One of the properties we expect 
from our quadrupole-aligned waveform is that during the in- 
spiral the frequency of the (I — 2. \m\ — 2) modes will to a 
good approximation satisfy the relation 



0)22 = 2(o) or b — cpcos 9) 



(1.1) 



Our main results are that (1) we can determine the 
quadrupole-aligned direction from the GW signal to high ac- 
curacy (within a fraction of a degree during most of the in- 
spiral), and (2) the GW signal is indeed far simplified, see in 
particular Fig. [TUJof the GW frequency before and after our 
(2,2)-maximization procedure, where the final frequency does 
indeed satisfy Eq. {OV In addition, we show that the GW 
signal is emitted in the direction of the orbital angular mo- 
mentum of the binary, which is not in general perpendicular 
to the orbital plane. We illustrate this effect with an example 
from PN theory, where it can be seen explicitly that the ef- 
fective orbital angular momentum is not parallel to the naive 
Newtonian angular momentum. 

In Sec. |ll] we provide details of our algorithm to find the 
orbital-angular-momentum direction from the GW signal, and 
in Sec.lIIIldescribe our numerical methods and numerical sim- 
ulations. The results of our method are presented in Sec. |IV] 
where we verify our method using a simple test case of an 
equal-mass nonspinning binary, and then apply the method 
to an unequal-mass spinning binary that undergoes significant 
precession. We discuss these results and prospects for future 
work in Sec. M 



II. MAXIMIZATION PROCEDURE ALGORITHM 

The Weyl scalar ^4 as calculated from the numerical code 
is decomposed into standard spin-weighted spherical harmon- 
ics (see [29] for our implementation). If the orbital angular 
momentum of the binary is parallel to the z-axis, then the GW 
signal will be dominated by the (I — 2, \m\ — 2) modes. We 
also expect that the coefficient of the (£ — 2, \m\ — 2) modes 
will be maximal in this case; for any other orientation of the 
orbital angular momentum, the (i = 2, \m\ = 2) modes will be 
weaker. 

Given the 1 = 2 modes W' 4 2m from the numerical code, we 
can rotate the frame to any other orientation using the transfor- 
mation described in Appendix [Aj to produce the correspond- 
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FIG. 1 : Profile of the magnitude of ^4,22 as the system is rotated by 
the Euler angles /3 and y. The example is taken from one time step 
(t = 562M) of the rotated equal-mass nonspinning case discussed in 
Sec. jlV A| Note that there is a clearly defined maximum, which in 
this case is at (j3, y) = (-10°, -205°). 



ing ^4.2™ in that frame. We locate the direction of the orbital 
angular momentum by searching over a range of the Euler an- 
gles (fi , y) to find a global maximum in ^4,22- 

The procedure in practice is as follows. We start our anal- 
ysis after the passage of the pulse of junk radiation. Since we 
extract the GW signal at either R ex = 90M or R ex = 100M, 
we take the start time to be at about t = 150M. At this time, 
we produce a first guess of the direction of L from the lo- 
cation of the black-hole punctures at that time. This pro- 
vides a guess (j3o,7o) of the Euler angles by which to rotate 
the system. Given this initial guess, we then search over a 
range of (j3,y) = (j3o ± 10°, yo ± 10°) with an angular res- 
olution of 0.1°, and find the angle for which the function 
l^^l 2 + | v I / 4,2— 2 1 2 has a maximum. In our test cases, where 
the orientation is constant, this procedure is trivial, but in gen- 
eral this first guess may not be very accurate. In particular, it 
does not take into account the time lag from the source to the 
GW extraction sphere. However, we do not expect the system 
to precess by as much as 10° over ss 100M of evolution. We 
also know that the Newtonian orbital angular momentum Ln 
calculated from the puncture motion is not in general parallel 
to the direction that maximizes the (£ = 2,m = 2) mode, but 
we do not expect the deviation to be larger than a few degrees; 



we will discuss this point further at the end of Sec. IV 



For subsequent times, we use the angles from the previ- 
ous time step as the first guess, and now search over the 
smaller range of ±3° in each angle. We locate the maximum 
in | v I / 4..22 1 2 + | v I / 4.2— 2 1 2 with a quadratic curve fit through the 
data from the search. 

At all times we find a clear maximum in the amplitude of 
¥4,22 as a function of the rotation angles. An example is given 
in Fig. [T] based on one time step of the rotated equal-mass 
nonspinning case presented in Sec. IV A 



III. NUMERICAL METHODS AND SIMULATIONS 

We performed numerical simulations with the BAM code 
[301 . The code starts with black-hole-binary puncture 
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initial data OTl [32l generated using a pseudo-spectral ellip- 
tic solver 11331 . and evolves them with the ^-variant of the 
moving-puncture |5j HO ED version of the BSSN ||35l l36l 
formulation of the 3+1 Einstein evolution equations. Spa- 
tial finite-difference derivatives are sixth-order accurate in the 
bulk QUI . Rreiss-Oliger dissipation terms converge at fifth or- 
der, and a fourth-order Runge-Kutta algorithm is used for the 
time evolution. The gravitational waves emitted by the binary 
are calculated from the Newman-Penrose scalar ^4, and the 
details of our implementation of this procedure are given in 
(29). See e.g. l37ll for a recent extensive parameter study of 
non-precessing binaries that uses the same numerical code and 
general setup. 

In each simulation, the black-hole punctures are initially a 
coordinate distance D apart, and are placed on the y-axis at 
yi = -qD/(l +q) and y 2 = D/(l +q), where q — M 2 /M\ is 
the ratio of the black hole masses in the binary, and we al- 
ways choose M\ < M 2 . The masses M; are estimated from 
the Arnowitt-Deser-Misner (ADM) mass at each puncture, 
according to the method described in |3"T1 ; see also the Ap- 
pendix of 11371 . The Bowen-York punctures are given mo- 
menta p x — TPt tangential to their separation vector, and 
p y = ±p r towards each other. The latter momentum com- 
ponent accounts for the (initially small) radial motion of the 
black holes as they spiral together. Initial parameters for low- 
eccentricity inspiral were produced using integrations of the 
PN equations of motion, as described in ll37l[38l . 

The eccentricity is measured with respect to the frequency 
of the orbital motion, as in all of our past work on eccen- 
tricity removal [ 37-40 1, and also discussed in IPTT1 l42l and 
references therein. The eccentricity is estimated as the ex- 
trema of e a (t) = (co(t) - C0Qc(t)) / (2o)Qc{t)), where (0 is the 
frequency of the (£ = 2,m = 2) mode of the waveform, and 
COQc(t) is an estimate of the frequency evolution for a non- 
eccentric binary, calculated by a smooth curve fit through the 
numerical data. 

The grid setup is similar to that used in [29|, and using the 
notation introduced there, the simulations discussed in this pa- 
per all use a configuration of the form XMr\=2 [h xiVifex 2N : 
6] . This indicates that the simulation used the % variant of the 
moving -puncture method, l\ nested mesh-refinement boxes 
with a base value of A/ 3 points surround each black hole, and 
I2 nested boxes with (2A/) 3 points surround the entire system, 
and there are six mesh-refinement buffer points. The 77 param- 
eter in the BSSN system is Mr\ = 2. The choices of N, l\, l 2 
and the resolutions are given in Tab. [I] The resolution around 
the puncture is denoted by M\jh m \ n , which is the resolution 
with respect to the smallest black hole, M\. The puncture of 
the second black hole will have the same numerical resolu- 
tion, but if the black hole is bigger, M 2 > M\ , then it will ef- 
fectively be better resolved. In unequal-mass cases, different 
numbers of refinement levels can be used around each black 
hole, so that the larger black hole need not be unnecessarily 
well-resolved, which would slow down the code. 

Far from the sources, the meaningful length scale is the total 
mass of the binary, M = M\ + M 2 , and so the resolution on the 
coarsest level is given by h mux /M. 

We consider two configurations. The first is an equal-mass 



nonspinning binary, using the same setup as first described 
in |43|. The initial separation is D = 12M, and the binary 
completes about nine orbits before merger. One additional 
simulation was performed in which the orbital plane was first 
rotated by 10° about the y-axis, and then around the z-axis by 
25°. 

The second configuration is a binary with mass ratio q = 3, 
where the larger black hole has spin 5 2 /M| = 0.75. In the 
calculation of the initial parameters, the spin is directed per- 
pendicular to the orbital angular momentum when the binary 
is at a separation of D = 30M. The configuration is evolved 
using the PN equations of motion until about D = 10M, and 
the momenta read off from the PN evolution at a point where 
the point particles pass through the xy plane. This leads to the 
parameters given in Tab. [I] 

Some key physical properties of the simulations are given 
in the last three rows of Tab.|I] the estimate of the eccentricity 
of the binary, the time when the GW signal reaches its peak 
amplitude, and the number of GW cycles up until that time. 

TABLE I: Parameters for the two configurations that we con- 
sider in this paper, the equal-mass nonspinning case, and the q = 3 
precessing-spin case. (For the rotated equal-mass nonspinning case, 
the momenta are p, = +{0.07567,0.03588,0.01477}.) The lower 
rows of the table indicate the numerical grid, which follows the con- 
ventions of (29] [37). 



q 


1 


3 




{0.488278,0.488278} 


{0.4779361,1.0234487} 


Si 


{0,0,0} 


{0,0,0} 


s 2 


{0,0,0} 


{-1.048,1.199,0.560} 


Xl 


{0,6,0} 


{0,15.0779,0} 


X2 


{0,-6,0} 


{0,-5.02598,0} 


D/M 


12.00 


10.05 


Vx 


^0.085035 


+"0.126292 


Py 


±0.000537 


+0.00139578 


Pz 





±0.0696932 


N 


64 


96 


(h,h) 


(5,5) 


(4/5,8) 


Mi/h mia 


21.3 


60.0 


hmzx/M 


12 


17.06 


Xi,mm/M 


774 


1647 


e 


0.0016 


0.0015 




1940 


1170 


cycles 


19 


14 



IV. NUMERICAL RESULTS 

A. Test case: equal-mass nonspinning binary 

In order to test our maximization procedure, we consider 
two simulations of an equal-mass nonspinning binary. In one, 
a reference case, the orbital angular momentum is oriented 
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FIG. 2: Motion of one black-hole puncture for the reference and 
rotated cases. The orbital planes are related by a rotation about the 
y-axis of 10°, and about the z-axis of 25°. 



the inspiral, for example, ^4 21 1 was reduced from ~ 10 4 to 

~io- 7 . 

These results demonstrate that our method works, and give 
us an indication of the error bounds. We expect that in gen- 
eral the errors will depend on the orientation angles of the 
binary, and will be worse when the angles are small. In these 
cases the sub-dominant modes will be smallest, and therefore 
will be resolved with less accuracy in the numerical code, and 
will then contribute more noise to the waveform in the rotated 
frame. However, we will take the errors from this example 
as the basis for our error bounds in other applications of our 
method. 



parallel to the z-axis, and so the waveform is already in the 
quadrupole-aligned frame. The simulation starts at D = 12M 
and covers about nine orbits before merger. 

In the second simulation the orbital plane is rotated. The 
orbital plane is first rotated about the y-axis by 10°, and then 
around the z-axis by 25°. The motion of the punctures in both 
the reference and rotated cases is shown in Fig. [2] The modes 
of *P4 g m are now mixed, and the power in the ^4,22 mode is 
distributed amongst the other (£ — 2) modes. This can be seen 
in Fig. bl In the reference case (denoted by 'i , 4/ m ), the (I = 
2,m = 1) mode is zero by symmetry, and the (jt = 2,m = 0) 
mode is dominated by numerical noise. In the rotated case, 
however, both sub-dominant modes have become significant. 
Note that oscillations are visible in the (£ = 2,m = 0) mode 
amplitude because it is a purely real function. 

We now want to see if our maximization procedure can 
recover the waveform from the reference simulation. In our 
procedure we search for a rotation of the system by the Euler 
angles (J3, y) such that the coefficient of the (£ = 2, \m\ = 2) 
modes is maximized. If the method works, we will recover the 
Euler angles (—10°, —205°), which correspond to the rotation 
we have described 1 . 

Fig |4] shows the error in the determination of the Eu- 
ler angles. The maximization procedure was applied from 
t = 150M, after the burst of junk radiation has passed, through 
to t = 2000M, which is late in the ringdown phase. Up until 
about t = 500M the waveform is rather noisy, and so the er- 
ror in /3 can be as large as 1°, and in y the error is up to 4°. 
During most of the inspiral, however, when the wave signal 
is clean, the error in /3 is below 0.05°, and the error in y is 
below 0.2°, and even during ringdown (when the waveform 
amplitude is falling exponentially), the angles are determined 
to within ±(0.5°, 2.0°). 

The magnitudes of the 1 — 2 modes in the quadrupole- 
aligned waveform agree well with those in the reference case. 
The (I = 2, \m\ — 2) modes agreed within numerical error in 
the raw data, and the {I = 2, \m\ — 1) modes, which should be 
zero by symmetry, were reduced by three orders of magnitude, 
to a level that would generally be regarded as noise. During 



The Euler angle to reverse the twist is —205° due to the freedom in per- 
forming the rotation about the y-axis clockwise or counterclockwise. 



B. Precessing binary 

Having shown that the maximization procedure works for 
the equal-mass nonspinning test case, where the orientation 
of the orbital plane is known, we now apply the method to 
a precessing binary. The configuration we have chosen has 
a mass ratio of q = M2/M1 = 3, the larger black hole has a 
spin of S2/M? = 0.75, and the spin initially lies in the orbital 
plane , i.e., perpendicular to the Newtonian orbital angular 
momentum. The small black hole is not spinning. 

We expect this configuration to exhibit significant preces- 
sion. The leading post-Newtonian contribution due to spin is 
the spin-orbit interaction, which can be characterized by the 
Hamiltonian 11441 (see e.g. also B31 ) 

H so =2^^, (4.1) 

where r is the coordinate separation of the black holes, and 
the effective spin S e ff, which is defined as 

where in our case one of the spins would be zero. From the 
spin-orbit interaction one can derive a post-Newtonian evolu- 
tion equation for the black-hole spin [44], 

S = -^S eff xL. (4.3) 

This indicates that the maximum amount of spin precession 
will be achieved when the spin is perpendicular to the orbital 
angular momentum. If one of the black holes has a Kerr pa- 
rameter Sj/Mf, then S will be largest if the larger black hole 
is spinning. This is also convenient from a numerical point 
of view, because the resolution requirements increase both as 
the mass is decreased, and spin is added; it is computationally 
cheaper to put the spin on the larger black hole. 

We also know from PN theory that S = — L to leading order. 
If we increase the mass ratio, then the orbital angular momen- 
tum L at a given separation will decrease, but the magnitude 
of the spin will stay the same. Therefore the relative change in 
L due to the precession of the spins will increase. This means 
that we will get greater spin precession for higher mass ratios. 
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FIG. 3: The left panel shows the amplitude of the ^n^,m modes for the reference case. The (£ = 2,m = 1) mode is zero by symmetry, and we 
see that the [i = 2,m = 0) mode is much smaller than the dominant mode, and is essentially noise during most of the inspiral. The right panel 
shows the corresponding amplitudes for the rotated case. We now see that both sub-dominant modes have become significant. The amplitude 
of the (£ = 2,m = 0) mode is oscillatory because it is a purely real function; see text for more details. 
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FIG. 4: Error in the angles for the tilt (/3) and twist (y) of the orbital 
plane, as determined by the maximization procedure. 



We have chosen q — 3 because this is reasonably large com- 
pared to typical simulations we have performed in the past, 
but low enough that we still expect to be able to achieve high 
accuracy. 

Fig. [5] shows the orbital motion of the two punctures in the 
simulation. The precession of the orbital plane is clearly visi- 
ble in the figure. 

Considering the leading order spin-orbit interaction 




FIG. 5: Motion of the black-hole punctures for the q = 3 precession 
simulation. The motion of the small black hole is shown in red, and 
the large black hole is shown in black. The precession of the orbital 
plane is clearly visible through the inspiral. 



naries. The time evolution of the momentum vector p is given 
by the Hamiltonian evolution equation 

dp _ dH 
dt dr 

If the Hamiltonian H depends on the spins, then consequently 
the momentum also picks up a contribution from the spins, 
and the velocity vector r is in general not parallel to the mo- 
mentum p. Consequently, the directions of the orbital fre- 
quency vector £2, 



(4.4) 



(4.5) 



Eq. (4.1 1 also exhibits another subtle feature of spinning bi- 



is in general not aligned with the angular momentum L = r x 
p. For the spin-orbit interaction defined by the Hamiltonian in 
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Eq. ( |4.1| i, this contribution to the angular momentum can be 
computed as B4l 



M ( / dm . 
-nx nx ( 3S+— A 



Ls °-M 



r x ( vx ( s+ ¥ A 



where 



(4.6) 



(4.7) 



and Lns i s tne nonspinning contribution to the angular mo- 
mentum (which is parallel to the vector rxv),v = r, and n is 
the unit vector in the direction of r. 

Note that the effect of non-alignment of £2 and L is maximal 
when the spin S is in the orbital plane. This is indeed the 
case for our initial conditions, and also during the evolution 
the spin component out of the orbital plane is significantly 
smaller than the components in the orbital plane. Note also 
that since the spin typically varies on a timescale larger than 



the orbital time scale, Eq. (4.6 1 will lead to oscillations in the 
angle between £2 and L with roughly the orbital period. 

Such oscillations are not present in the direction of L, as 
illustrated in Fig. [6] We will see later that the quadrupole- 
aligned frame moves consistently with L (i.e,. as a smooth 
function), suggesting that our maximization procedure tracks 
the direction of the orbital angular momentum. 

The left panel of Fig. [7] shows the amplitude of the (t — 
2,m = 2) and (£ = 2,m = 1) modes during the inspiral. We 
clearly see that the "sub-dominant" (2,1) mode is of com- 
parable magnitude to the (2,2) mode, and shows significant 
modulation. (It is also instructive to compare with the re- 
sults in J46), where a precessing binary is also considered 
from a fixed frame of reference, and all of the £ = 2 modes 
are of significant amplitude.) The right panel of Fig. [7J shows 
the frequency of the (2,2) mode, ©22 = 022, over the same 
time interval. The frequency clearly exhibits large oscilla- 
tions. Based on the discussion around Eq. (jTTTJ we expect 
oscillations in ©22 of purely physical origin, but we also as- 
sume that the physical oscillations will be exaggerated and 



their frequency modified in the fixed frame of an inertial ob- 
server. 

We now apply the maximization procedure to the waveform 
signal from t = 150M, when the junk radiation has passed, 
through merger and ringdown (up to t = 1250M). At each 
time step the system is rotated such that the {I = 2, \m\ =2) 
mode amplitudes are maximized. 

Finally, we address the question of whether the GW signal 
is emitted normal to the orbital plane, or parallel to the orbital 
angular momentum. Although we cannot unambiguously de- 
fine the direction of orbital angular momentum, we can cer- 
tainly determine whether the GW signal is emitted normal to 
the orbital plane. 

Fig. [8] shows the Euler angles (j3, 7) that were found in the 
maximization procedure, time shifted by 100M to approxi- 
mately compensate for the time lag to the extraction spheres. 
It also shows the angles (9,(p) of the direction orthogonal to 
the orbital plane as computed from the NR simulation, and for 
the orbital angular momentum L as computed from a PN sim- 
ulation (as in Figs. [6]). A constant offset of 2° has been added 
to 9 to align the 9 angles and j3 . A potential explanation for 
this small offset are coordinate gauge ambiguities. If the GW 
signal were emitted normal to the orbital plane, we would ex- 
pect to be able to align j3 with —9 from the numerical relativ- 
ity simulation, and likewise for 7 and — (p. However, it is clear 
from Fig. [8] that the orbital-plane angles contain extra oscil- 
lations. Based on the illustration in Fig. |6| this suggests that 
the GW signal is emitted in the direction of the orbital angular 
momentum. In particular, show in Fig. [8] the direction of the 
orbital angular momentum as predicted in PN theory shows 
good agreement with the angles that define the quadrupole- 
aligned frame. 

Fig. |9j shows the amplitude of the original ¥4 2 2 an d me 
quadrupole-aligned signal that results from the maximization 
procedure, ¥4 22- We see that the maximization procedure has 
indeed increased the amplitude at all times, and also seems to 
have removed some oscillations. 

The frequency of the (£ = 2,m = 2) mode before and after 
the maximization procedure is shown in Fig. 10 This figure 



illustrates the key result of this work: the high-frequency os- 
cillations in the wave frequency have been removed by the 
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FIG. 8: The Euler angles (/3 , y) found when the maximization procedure was applied to the q = 3 precessing-binary waveform. For comparison 
we show the corresponding angles (— 6,— <p) of the normal to the orbital plane as computed from the NR simulation, and for the angular 
momentum L from a PN simulation (as in Figs. [6}. We apply an appropriate time shift to (/3,y) to approximately compensate for the time 
lag at the wave extraction sphere, and add 2° to the curves for as described in the text. We clearly see that the orbital-plane angles show 
additional oscillations that are not present in the (2,2)-maximization angles. 



maximization procedure, and we are left with a far simpler 
functional form. We note, however, that the oscillations in the 
frequency have not been completely removed. This is to be ex- 
pected from Eq. < | 1 - 1 J > . In the absence of precession, during the 
inspiral the gravitational wave frequency of a spherical har- 
monic mode (£, m) is with a high degree of accuracy propor- 
tional to the orbital frequency, (0( m — m(0 ol ±,. In the presence 
of precession, this is however replaced by Eq. (JTTTJ, which 
adds an extra term depending on the precessing motion of the 
orbital plane. In Fig. 11 we compare the frequency of the 
(£ = 2,m = 2) mode after the maximization procedure with 
the orbital frequency with the precession term added accord- 
ing to Eq. and we find reasonable agreement. We also 
show the frequency co^ that results from rotating the system 
according to the direction perpendicular to the orbital plane, 
which is also the direction of the naive Newtonian orbital an- 



gular momentum. It is clear from Fig. 1 1 that the oscillations 
due to the orbital-plane rotations are much larger, and this fur- 
ther suggests that the quadrupole-aligned frame is optimal. 

It is clear that the maximization procedure produces (£ — 
2, \m\ = 2) modes that are of a simpler form than in the origi- 



nal data. However, this is not a guarantee that we have cor- 
rectly tracked the direction of the GW emission; we have 
not necessarily put the waveform into a physically meaning- 
ful frame of reference. One test of our method is to calcu- 
late the effect on the sub-dominant modes. We expect that in 
the quadrupole-aligned frame the amplitude of the GW sig- 
nal will agree to a good approximation with that from a q = 3 
nonspinning binary. (The spin effect on the rate of inspiral is 
dominated by S • L, and this is close to zero throughout our 
simulation, so we expect the inspiral to be similar to that for a 
nonspinning binary with the same mass ratio.) 



Fig. 12 shows a selection of modes for the quadrupole- 
aligned waveform. The left frame shows the transformed 
modes for the precessing binary, and the right frame shows 
the same modes for the nonspinning q = 3 waveform pre- 
sented in 11371 . Two things are remarkable about this figure. 
The first is that the amplitudes of the modes show extremely 
good agreement. The other is that we have found that the 
magnitude of the (£ = 2,m = 1) mode is extremely sensitive 
to the angle by which the system is rotated. If, for example, 
we were to modify j5 or y by a fraction of a degree, ^4 21 
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600 700 800 900 1000 1100 

t[M] 

FIG. 10: Frequency of the (£ = 2,m = 2) mode before QSf^ 2 2> ar, d 
after (¥422) me maximization procedure. We see that the high- 
frequency oscillations have been removed. The remaining oscilla- 
tions are of a lower frequency and much lower amplitude; see text 
and Fig. 1 1 1 1 

could change by orders of magnitude. With this fact borne in 
mind, the oscillations in ^4,21 | are not very large at all. This 
figure suggests that we have located an optimal frame from 
which to study the GW signal. 



V. DISCUSSION 

We have presented a simple method to track the precession 
of a binary system, using only information from the GW sig- 
nal. Our procedure is to rotate the system such that the mag- 
nitude of the {I = 2, \m | = 2) modes is maximized, based 
on the physical assumption that this will be the direction of 
dominant GW emission. We refer to this as the "quadrupole- 
aligned" waveform. Based on evidence from PN theory, we 
show that this direction seems to correspond to that of the or- 
bital angular momentum, which is in general not perpendic- 
ular to the orbital plane. We also show that our method pro- 
duces higher-mode amplitudes consistent with what we know 
from non-spinning, non-precessing binaries. 

The result of our procedure is that the waveform is repre- 
sented in a more simple form than the one produced directly 



500 600 700 800 900 1000 

t[M] 

FIG. 11: Frequency of the (I = 2,m = 2) mode after <^Ha,22) the 
maximization procedure, compared with the orbital frequency with 
a precession term added according to Eq. l |l.l[ >. We also show the 
frequency that results from rotating the system according to the di- 
rection of the Newtonian orbital angular momentum, (%, i.e., the 
normal to the orbital plane. 



from the numerical code. This is particularly true of the sub- 
dominant modes; compare Figs.[7]and[T2] We expect that this 
will simplify the task of producing analytic inspiral-merger- 
ringdown models, which is one of the main motivations for 
our work. This method also provides a normal form for the 
waveform, which will facilitate future comparisons between 
numerical and analytic results. 

One could propose alternative procedures to track the pre- 
cession of the system, and we will now discuss some of them, 
and their difficulties. 

Only the total angular momentum of the spacetime is unam- 
biguously defined in General Relativity. The form of Bowen- 
York puncture initial data are such that we can analytically 
calculate the angular momentum ([29, 47 H491 ) of the initial 
slice from the initial-data parameters; it is simply given by 
L = ri x pi +T2 x P2, where r,- are the coordinate locations of 
the punctures, and p, are the momenta that are input into the 
Bowen-York extrinsic curvature. We can calculate the angular 
momentum radiated through the spheres on which we measure 
the GW signal, and so we can determine the total angular mo- 
mentum of the system as a function of time. However, we 
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FIG. 12: Left: selected modes of the precessing-binary waveform, after being transformed into the non-precessing frame, i.e., after the system 
has been rotated by the angles that were found from the (2,2)-maximization procedure. The right-hand plot shows the same modes for a 
nonspinning (and therefore non-precessing) q = 3 waveform. The agreement is remarkable. Note in particular the qualitative agreement of the 
(£ = 2,m = 1) mode, which is of comparable magnitude to the (£ = 2,m = 2) mode in the raw data (see Fig.[7]l. 



want the orbital angular momentum, L = J S. To calculate 
this we need to know the black-hole spins as a function of 
time (which can be estimated with reasonable accuracy from 
the black holes' apparent horizons [50]), but these quantities 
are calculated at the black holes, not at the GW extraction 
sphere, and cannot easily be translated. 

One could attempt to instead calculate the orbital angular 
momentum entirely at the sources, but this also presents diffi- 
culties. The proper distance between the black-hole horizons 
and their momenta could be calculated by some quasi-local 
procedure (for example BP ), and hence the orbital angular 
momentum. But it will be difficult to assess the gauge errors 
in any such method. Alternatively, one could calculate the an- 
gular momentum using the puncture locations and PN theory, 
but this will only be an approximation to the true general rel- 
ativistic angular momentum. One direction we can easily de- 
termine is the normal to the orbital plane of the binary, but we 
have seen in Sec. |IV| that this is not the direction in which the 
dominant GW signal is emitted, and nor does it define a ref- 
erence frame from which the GW signal appears simpler than 
what can be achieved by the maximization procedure that we 
have used. 

Nonetheless, a number of issues remain to be resolved 
in our procedure. In particular, our method does not seem 
to accurately track the quadrupole-aligned direction through 
merger and ringdown. If it were able to do this, it would pro- 
vide an alternative procedure to determine the direction of the 
spin of the final black hole. We find that the angles from our 
maximization procedure continue to vary through merger and 
ringdown, and do not settle at constant values, which is what 
they would do if they had determined the final spin direction. 
This problem may be due to the accuracy of the numerical 
data, or to more subtle effects, for example the motion of the 
center-of-mass of the system due to gravitational recoil. We 
will investigate this issue further in future work. 

While working on this project we have learned that an in- 
dependent effort to identify precession effects via a similar 
algorithm will be presented by Seiler et al. in a forthcoming 
publication (52] • 
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Appendix A: Transformation of ^am under rotations 

We aim to derive the transformation of the Weyl scalar ¥4 
under a rotation R 6 SO (3). It can be shown that the Weyl 
scalar is a field of spin-weight s = —2 and hence it can be 
expanded as 

^4 = £^4,/mI'te , (Al) 
l.m 

where F /m 2 denote the spherical harmonics of spin-weight s — 
—2 1531 . For s = we obtain the regular spherical harmonics 
Y[ m , which are the eigenfunctions of the angle-dependent part 
of the Laplace operator. 

The transformation of the spin-weighted spherical harmon- 
ics is a simple composition of the transformation of the spin- 
basis-dependent part and of F/ m . It is convenient to introduce 
standard polar coordinates (r,0,(p) and to define Y\ m with re- 
spect to the polar angles (9, <p). The spherical harmonics then 
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have the form 



the following law under a change of the spin basis: 



(A2) 



We will consider rotations R, which transform angles £2 = 
(9,(p) to the new coordinates £2' = (0', <p'). The spin-weight- 
zero spherical harmonics l/ m then transform according to 
Yi m (8,(p) i — y Yi m (Q',(p') by applying the operator Pr, where 
R is a rotation about the z-axis by the angle y such that 
<p h->- <p' = <p + 7 and = 0', is given by 

Y lm (e',(p') =p R Y lm (e,(p) = ^Y ln {e,(p). (A3) 

Now, let R(y/3a) denote an arbitrary rotation by the Euler 
angles 7,j3,o:. Using the z-y-z convention, the spherical har- 
monics then obey the following transformation law ll54l 1551 : 



Y lm (8',<p')= £ ^dl^y^YUO,^ (A4) 

m'=-l 

where the d 1 , denote the Wigner of-matrices which are given 

mm fe fe 

by 

d' m , m =^(l + m)!(l-m)\(l + m>)l(l-m>)\ 

^ ^ ^k-\-m'—m 



Y k\{l+m-k)\(l-m> -k)\(m> -m+k)\ 

X (Sin ^) 2k+m '-' n {cOS ^f-2k-m'+m_ 



(A5) 



Due to the properties of the group SO(3), the inverse transfor- 
mation is then given by 

Y lm (9,(p)= £ e- im 'rd l m , m (-P)e- ima Y lm ,(9',cp'). (A6) 

m'=-l 

The next step is to include the change of spin-basis under a ro- 
tation. According to [56 1 a quantity rj of spin-weight s obeys 



77 =77 



(A7) 



Combining Eqs. ( A6 1 and ( A7 1 yields the transformation law 



for the spin-weighted spherical harmonics: 



Yf m (6,cp)=e-^ £ e-'-rd'^J-^e-^Y^e'^'). 

m'=-l 

(A8) 



We invert Eq. ( Al 1 to determine the transformation law for 
the y i l 4j m -modes, 



x t J 4j m = I ¥4Y? m (e,(p)dQ. 



xe" na Y l s m ,(9',( P ')dQ.' 

m'=-l 

where we see that explicit knowledge of J as a function of 9 
and <p is not necessary to determine the coefficients l i , 4/ m . 
This transformation law can now be applied to any given 
lm , e.g., our numerical data, in order to change the frame of 
reference. The remaining free parameters are the three angles 
that determine the rotation. We restrict ourselves to a rotation 
about two Euler angles, j3 and 7, only. Since we aim to align 
the orbital angular momentum with the z-axis at every instant 
of time, i.e., L 1— > z, a simple calculation shows that in order to 
fulfill this j3 = — 9 and 7= — <p are required, where (9, (p) are 
the polar coordinates determining the direction of L. 
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